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ABSTRACT 

Millisecond and binary pulsars are the most stable natural frequency standards which 
allow the introduction of modified versions of universal and ephemeris time scales 
based on the intrinsic rotation of pulsar and on its orbital motion around the barycen- 
tre of a binary system. The measured stability of these time scales depends on numer- 
ous physical phenomena which affect the rotational and orbital motion of the pulsar 
and observer on the Earth, perturb the propagation of electromagnetic pulses from 
the pulsar to the observer, and bring about random fluctuations in the rate of time 
measured by an atomic clock used as a primary time reference in timing observations. 
Over long time intervals the main reason for the instability of the pulsar time scales 
is the presence of correlated, low-frequency timing noise in residuals of times of ar- 
rivals (TOA) of electromagnetic signals from the pulsar which has both astrophysical 
and geophysical origins. Hence, the timing noise can carry important physical infor- 
mation about the interstellar medium, the interior structure of the pulsar, stochastic 
gravitational waves coming from the early universe and/or binary stars, etc. Each 
specific type of low-frequency noise can be described in a framework of the power law 
spectrum model. Although data processing of pulsar timing observations in the time 
domain seems to be the most imformative, it is also important to know to which spec- 
tral bands single and binary pulsars, considered as detectors of the low- frequency noise 
signal, are the most sensitive. A solution to this problem may be reached only if the 
parallel processing of timing data in the frequency domain is fulfilled. This requires a 
development of the Fourier analysis technique for an adequate interpretation of data 
contaminated by the correlated noise with the noise spectrum which is divergent at 
low frequencies. The given problem is examined in the present article. 
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1 INTRODUCTION 



It is well known that millisecond and binary pulsars provide an exellent means for testing the 
theory of general relativity (Taylor & Weisberg 1982, 1989, van Straten et al. 2001, Weisberg & 
Taylor 2003, Kramer et al. 2003). Pulsar timing is a powerful tool for exploring the structure of 
the interstellar medium (Rickett 1990, 1996, Shishov 2002) and for investigating the interior of 
neutron stars (Cordes & Greenstein 1981, Kaspi et al. 1994) as well as setting an upper limit 
on the energy density of gravitational radiation produced in the early universe and/or by the 
orbital motion of binaries (Kaspi , Thorsett & Dewey 1996, McHugh et al. 1997, Kopeikin 1997a, 
Kopeikin 1999a, Kopeikin & Potapov 2000, Lommen 2002, Jaffe and Backer 2003). Physical process 
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of pulsar's intrinsic rotation around its axis has been proposed and used (Backer et al. 1982, Il'in 
et al. 1986, Rawley et al. 1987, Guinot & Petit 1991, Matsakis et al. 1997) as a new time reference 
(PT scale), analogous to the universal time (UT) in time-keeping metrology. The UT is based on 
observation of the Earth's diurnal rotation. The rate of the atomic time scale UTC is regulated 
by introduction of leap seconds to reproduce the rate of UT and eliminate the secular difference 
between the uniform atomic time scale TAX and UT (Kovalevsky and Seidelmann 2004). Quite 
recently, a new astronomical time reference (BPT scale) spread out over extremely long intervals 
(approaching 50-100 years) has been suggested (Petit & Tavclla 1996) and its stability have been 
carefully examined by Ilyasov et al. (1998) and by Kopeikin (1999). BPT scale is derived from the 
orbital motion of a pulsar in a binary system and represents an analogue of ephemeris time (ET) 
in classical astronomy as based on the orbital motion of the Earth around the Sun (or the Moon 
around the Earth) and introduced by Newcomb (1898) at the end of the last century. 

An adequate analysis of the timing data requires a deeper understanding of the nature of the 
noise process present in pulsar timing residuals. As soon as the autocovariance function of the noise 
process is known, analysis in the time domain becomes possible. Time domain analysis is the most 
informative since the observed stochastic process is not stationary and includes a non-stationary 
component which contribution to the parameters' estimates can be important (Groth 1975, Cordes 
1978, 1980; Kopeikin 1997b). Because pulsar timing observations arc conducted over relatively long 
time intervals, the white noise from errors in measuring TOA of pulsar's pulses is suppressed by the 
presence of a number of correlated, low-frequency ( "red" ) noise processes having different spectra 
and intensities. Henceforth, we are mainly interested in analysing the red noise. 

A simple, but realistic mathematical model for the red noise has been developed by Kopeikin 
(1997b). The model is based on the shot noise approximation with autocovariance function de- 
pending on both stationary and non- stationary components. The given model elucidates a rather 
remarkable fact (Kopeikin 1999b), namely, that the timing residuals as well as variances of some 
spin-down and all orbital parameters are not affected by the non-stationary component of the red 
noise at all because its contribution is filtered out by the fitting procedure based on the least square 
method. This observation justifies and puts on a firm ground the Fourier analysis of TOA residuals 
and variances of fitting parameters in the frequency domain which is normally done under implicit 
assumption that the non-stationary component of the autocovariance function is either absent or 
unimportant. We conclude that the Fourier analysis gives unbiased information about the noise 
process itself and permits us to reveal which frequency harmonics in the spectral expansion of the 
stochastic process of the pulsar timing observations are the most substantial. 

Any low-frequency noise can be approximately characterized by the power-law spectrum S{f) ~ 
where the spectral (integer) index n ^ 1. It is obvious that the spectrum has a singularity at 
zero frequency. Hence, the energy of TOA residuals comprised in such noise should formally have 
infinite value because the integral over all frequencies from zero to infinity taken from S{f) ~ 
is divergent. Clearly, this has no physical meaning, and one has to resort to special mathematical 
methods in order to avoid this illicit divergency. There are two standard mathematical methods 
for curing this flaw. The first method consists in the analytical continuation of the spectrum by 
formally changing it form from S{f) ~ to S{f, A) ~ j^-(n+A) -^j^gj-g A is a purely imaginary 
parameter different from zero. Such a model of the analitically continued spectrum gives convergent 
integrals which coincide everywhere on the real axis with the integrals from the original spectrum 
S{f) ~ /""", exept for the point / = 0. In order to prescribe a physical meaning to such integrals, 
we have to expand them in a Laurent series with respect to the parameter A and take the finite 
part of the expansion for ^4 = 0. Such procedure has been used, in particular, by Kopeikin (1997a) 
for regularization of the autocovariance function of stochastic noise of the primordial gravitational 
wave background having formally divergent spectrum S{f) ~ (Sazhin 1978, Detweiler 1979, 
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Mashhoon 1982, 1985, Bertotti et al. 1983, Mashhoon & Seitz 1991) and for setting an upper limit 
on the energy density of this background noise in the ultra-low frequency domain. 

The procedure of analytic continuation of divergent integrals is mathematically rigorous and 
represents theoretically a powerful tool which gives well-defined and self-consistent mathematical 
results (Gel'fand & Shilov 1964). However, for researchers doing numerical computations the second 
method of regularization of singular spectra is preferable and more useful practically. It is based on 
the infrared cutoff of all divergent integrals at the frequency /cutoflf — where £ <S 1 is constant. The 
infrared cutoff technique requires corresponding modification of the power-law spectrum of the red 
noise in order to eliminate the dependence of the fitting parameters and the autocovariance function 
of the noise on the ad hoc introduced frequency cutoff e. We shall prove that this modification 
consists of subtracting from the original power-law spectrum the infinite series of the frequency- 
dependent Dirac's delta function and its derivatives with numerical coefficients which are functions 
of the cutoff frequency e. Such modification of the spectrum preserves the structure of the most, 
generally accepted, autocovariance functions and, as a consequence, does not change results of the 
numerical processing of signals in time domain. In addition, the infinite delta-function dependent 
part of the modified spectrum is rapidly convergent and accounting for contribution of few first 
terms is sufficient in practice. This second method of regulariazation of red spectra being divergent 
in the infrared frequency domain is worked out in the present paper. 

In what follows, it is more convenient for mathematical purposes to express results in terms of 
units in which frequency and time are dimensionless. For instance, in binary pulsars, it is preferable 
for symbolic calculations to measure time in units of orbital frequency, = 271 /Pi,, where Ph is 
the orbital period of the binary system. In these units, frequency / is measured in terms of 1/Pb, 
and dimensionless time is the pulsar's mean orbital anomaly u = ni,t where t is time measured in 
SI units. 

2 REGULARIZED SPECTRUM OF LOW-FREQUENCY NOISE 

Any gaussian low-frequency noise is completely characterized by the autocovariance function which 
describes the correlation between two values of stochastic process separated by the arbitrary time 
interval u = nh{t2 — ti). The autocovariance function usually consists of two algebraically indepen- 
dent parts characterizing separately a stationary, R~{u), and a non-stationary, /?"*"(«), component 
of the noise. Complete expressions of the autocovariance functions for various examples of low- 
frequency noise have been derived by Kopeikin (1997b), where the shot noise approximation of 
a stochastic process has been used. Although both stationary and non-stationary components of 
the autocovariance function are important for a comprehensive treatment of observations, we are 
dealing in the present paper only with the stationary part since the non-stationary component of 
the red noise is filtered out by the least square method (Kopeikin 1999b). 

The stationary part of the authocovariance function can be transformed into the spectral density 
of noise, S{f), by means of the Wiener-Khintchine theorem 



where u is the dimensionless time and / is the dimensionless Fourier frequency. If the (constant) 
intensity of noise is denoted by hn then the autocovariance function of low-frequency (correlated) 
noise is determined by the expression (Kopeikin 1997b) 




(1) 




CnhnU^ ^ In n — 1, 3, 5, fiicker noise 



n = 2, 4, 6, random walk noise 



(2) 
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where C„ is a numerical constant of normalization. 

Functions, which might be appropriate candidates for the spectrum of noise procesess with the 
foregoing autocovariance functions, are S{f) ~ hnf^^ where n is integer. However, integrals 
from such power-law functions are divergent because of non-physical singularity of the spectrum at 
zero frequency. For this reason, a regularization technique should be used because we don't usually 
know the asymptotic low-frequency behaviour of the spectrum and can not treat the extremaly 
low-frequency domain adequately. 

To develop the regularization method the power spectrum S{f) of the noise should be defined 
using the truncated Fourier transform of the autocovariance function with the lower cutoff frequency 
/cutoff = £■ Specifically, we require that the truncated cosine Fourier transform of S{f) must give the 
stationary part of the original autocovariance functions (0) without any additional contributions 
which might come up from the cutoff-frequency lower limit of the integral and bias the estimates of 
the fitting parameters. Let us guess that this requirement can be satisfied with the specific model 
of the red-noise spectrum S{f) represented by the formula 



Sif) 



1 oo 



0, 



k=0 



otherwise 



(3) 



where the spectral index of noise n = l,2,...,6, the constant parameter hn is the strength of noise, 
quantities Bk{e) are constant numerical coefficients being defined later, and 5^^\f — e) denotes the 
n — th derivative with respect to / of the Dirac delta-function 6{f — e). The Dirac delta function 
is defined according to the condition (Korn & Korn 1968) 

0, if X < a, or X > b, 



f{x)6{x — X)dx 



i/(^ + o), 
i/(x-o). 



if X = a, 
if X = 6, 



(4) 



^ i[/(X-0) + /(X + 0)], ifa<X<6, 



where f{x) is an arbitrary function being such that the unilateral limits /(X — 0) and /(X + 0) 
exist. The proposed form of the low-frequency spectrum S{f) will be acceptable if the coefficients 
Bk[e) can be uniquely determined by the condition that the cosine Fourier transform of S{f) gives 
the stationary part of the autocovariance function of a corresponding low-frequency noise shown in 
Eq. (j2)). This statement is valid and its consistency can be proven by making use of straightforward 
mathematical calculations. As an example, we show how to determine the several first coefficients 
Bk{e) in the event of fiicker noise in pulsar's rotational phase which has the spectral index n = 1. 
Coefficients Bk{e) for noises having other spectral indices will be displayed in this section without 
proof (which is rather easy). 

The stationary part of the autocovariance function R^{t) of the fiicker noise in the pulsar's 
phase is equal to hi7i~^ log \u\. According to the definition of the spectrum, accepted in the present 
paper, we should have 



R- 



2 Sif) cosi27i fu)df, 



(5) 



where the integral has been truncated at the cutoff frequency £ <^ 1 to eliminate its divergency 
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at the point / = 0. Substituting S{f) from Eq. Q with n = 1 on right hand side of Eq. (0) and 
taking the integral we get 

h °° 

R-{u) = iCi(27r£|M|) + hi COs(27r£M) V(-l)'=E2fc(£)(27r£M)2^ (6) 

^ fc=0 

where Ci(x) is the cosine integral, and we have used the formula (Korn & Korn 1968) 

poo /^2A; 

2 / 5{2fc)(j _ ^) cosi2n fu)df = —r cos(27reM) = {-lf{2nuf'' cos(27r£M). (7) 



A Taylor expansion of the cosine integral and cos(27r£M) in the right hand side of the expression 
(jH)) with respect to small parameter e yields 



R-(u) = -^1- log I u| — 7 — log(27r£:)] + hil Bo{e) + {2neu 



TT 



1 1 



An 2 

Since by definition R~{u) must be equal to —hiTi~^ log we find from Eq. 



TV 



7 log(27r£:) 



7r 



^ ^ Att 2tt 2n ' 



BJe)=0(e'). 



(9) 



We believe that for practical purposes it is enough to account for the coefficient -Bo(^) only, since 
all other residual terms that appear are multiplied by the factor 2neu, which is expected to be 
negligibly small under usual circumstances, because of the smallness of the product en. From the 
formally mathematical point of view the smaller the cutoff frequency e the more exact is our 
approximation. However, one should take into account that the real spectrum of the low-frequency 
noise can have the low-frequency behaviour different from the power law and our approximation 
will not get better for a very small value of e. In other words the residual terms in the model of the 
spectrum under discussion are model dependent. Had we chosen another model for the spectrum 
having different behaviour as the frequency approaches zero the residual terms would look different. 
What model of the low frequency noise works better in the procedure of fitting parameters of the 
pulsar is an interesting question for further discussion. We emphasize however that the definition 
of the low-frequency noise spectrum in terms of the power-law model amended by the series of 
delta-functions and its derivatives makes the estimates of the fitting parameters not sensitive to 
the cutoff frequency and, thus, extrapolates the power-law model of the noise spectrum to lower 
frequencies than the power-law model of noise which does not account for the delta functions. 

Proceeding in the same way for the set of other spectral indeces we obtain the following expres- 
sions for the power spectra of low-frequency noises: 

(1) Flicker noise in phase: 

SU) = ^{7 + 2 [7 + log(2vre)] 6{f - .)} + 0{e\ (10) 

(2) Random walk in phase: 

^(/) = ^{4-'^(/-^)}+o(^)- (11) 
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(3) Flicker noise in frequency: 



S{f) 



{73 - ^^(f - ^) + [log(2^^) + ^ - 1] '^^'H/ - e)] + 0{e'). 



(4) Random walk in frequency 



(12) 



S{f) 



hi 



1 



0{e). 



(13) 



(5) Flicker noise in frequency derivative: 



1 

12 



7 + log(27r£:) - 



5W(/-e)}+0(52).(14) 



(6) Random walk in frequency derivative: 



15^3 



5(2) (/- 



0(.). 



(15) 



Note that the coefficient Bi[e) = in the expression (fT3j) for the spectrum of the random walk in 
the pulsar's frequency derivative. 

The expressions given above indicate that there is a strong concentration of the infinite energy 
of the noise at the lower cutoff frequency. As we have stressed already it is a specific feature of 
the chosen model of the spectrum which appears because we do not know the real behaviour of 
the spectrum while the Fourier frequency is approaching zero. Another remark is that a formal 
integration of any of the foregoing spectra with spectral indices n ^ 2 with respect to frequency 
from f = e to infinity (which may be naively treated as the energy being stored in TOA residuals) 
gives zero value which may look surprising ^. However, it is worth noting that the whole amount of 
energy presents in TOA residuals can be calculated only after multiplication of the spectrum by the 
filter function of the fitting procedure (see section 6 below for more detail). Therefore, calculation 
of the total energy of residuals is more complicated and always gives a positive numerical value 
as expected. Similar arguments can be used in calculating variances of the fitting parameters. 
For example, the calculation of variances of the ffist several spin-down parameters in the Fourier 
frequency domain may give a negative numerical value for the variance (Kopeikin 1999b) which is 
physically meaningless. The paradox is solved if we take into account the contribution of the non- 
stationary part of the noise which always makes variances of the parameters numerically positive 
2 

Pulsar timing observations can be used for the estimation of the strength and spectrum of low 
frequency noise present in TOA residuals. For this reason, the development of practically useful 
estimators of spectrum of noise are required. We are not going to consider in the present paper 
the question of how to construct the best possible estimators. This subject has been significantly 
explored by a number of other researches (see, for instance, the papers of Deeter & Boynton (1982), 
Deeter(1984), Taylor 1991, Matsakis al. 1997, Scott et al. 2003). Our purpose is to study the spectral 
dependence of TOA residuals and the variances of fitting parameters which are used in real practice. 



^ A simplest way to see this result is to notice that such integration corresponds to the value of the autocovariance function R(u) taken 
at the point u = 0. 

^ For more detail see (Kopeikin 1999b) 
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In order to put our approach on a systematic basis connected to our previous works let us describe, 
first of all, the timing model we are dealing with. 



3 TIMING MODEL 

We consider a simplified, but still realistic model of arrival time measurements of pulses from a 
pulsar in a binary system. It is assumed that the orbit is circular, and the pulsar rotates around 
its own axis with angular frequency Up which slows down due to the electromagnetic (and other) 
energy losses. It also takes into account that the orbital frequency of the binary system, Uf,, and 
its projected semimajor axis, x, have a secular drift caused by radial acceleration of the binary 
(Damour & Taylor 1991, Bell & Bailes 1996), its proper motion in the sky (Kopeikin 1996), and 
the emission of gravitational waves from the binary (Peters & Mathews 1963, Peters 1964) bringing 
about the gravitational radiation reaction force (Damour 1983a, Grishchuk & Kopeikin 1983). The 
timing model, which we work with, is described in full detail in our previous paper (Kopeikin 1999b) 
which particular notations we use in the present paper. Let us summarize those equations which 
will be useful for our anlysis. 

We assume that all observations of the binary pulsar are of a similar quality and weight. Then 
one defines the timing residuals r{t) as a difference between the observed number of the pulse, 
A/'"'''^, and the number J\f(t,9), predicted on the ground of our best guess to the prior unknown 
parameters of timing model, divided by the pulsar's rotational frequency u, that is 

(16) 

where 6 = {6a, a = 1,2, ...k} denotes a set of k measured parameters ^ which are shown in Tabled 
It is worth noting that hereafter we use for the reason of convinuence the time argument u = nf,t, 
that is the current orbital phase. If a numerical value of the parameter 6a coincides with its true 
physical value 6a, then the set of residuals would represent a physically meaningful noise e(t), i.e. 

r{t,6) = e{t). (17) 

In practice, however, the true values of parameters are not attainable and we actually deal with 
their least square estimates Therefore, the observed residuals are fitted to the expression which 
is a linear function of corrections to the estimates ^* of a priori unknown true values of parameters 
6a- From a Taylor expansion of the timing model we work with (Kopeikin 1999b), and the fact that 
r{t, 6) = e(t) one obtains 

r{t, 6*) = e{t) - YlMait, 0*) + 0{Pl), (18) 

a=l 

where the quantities Pa = = ^a~^a are the corrections to the unknown true values of parameters, 
and ipa{t,6*) = {dAf /d6a)Q^Q, are basic fitting functions of the timing model. 

In the following discussion it is convenient to regard the increments (3a as new parameters whose 
values are to be determined from the fitting procedure. The parameters (3a and fitting functions are 
summarized in Tabled with asterisks omitted and time t is replaced in accordance to our notations 
and conventions by the function u = ribt which is the value of the orbital phase taken at time t. We 
restrict the model to 14 parameters since in practice only the first several parameters of the model 



^ The numebr fc = 14 in the particular timing model we are dealing with (see the paper (Kopeikin 1999b) for fuller detail. 
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Table 1. List of the basic functions and parameters used in the fitting procedure. Spin parameters 
SNo jSujS v,5 v,5 V ,5 V , fit rotational motion of the pulsar around its own axis. Keplerian parame- 
ters 5x, (5(7, Sni fit the Keplerian orbital motion of the pulsar about the barycentre of the binary system. 
Post-Keplerian parameters 5 x,S x,5 x,5 nij,S Uf, fit small observable deviations of the pulsar's orbit from 
the Keplerian motion caused by the effects of General Relativity, radial acceleartion, and proper motion of 
the barycentre of the binary system with respect to observer. 



Parameter 



Fitting Function 



5Uo 



(92 = 



1 Sf 



1 5_u_ 

o 



Py — —5x sin (T — SfTx cos it 
Bg — -—Sx cos a + Scrx sin a 



I3q = (^ — 6x cos a + 5ni,x sin ct^ 

/3iQ = {^ — ^ ^ sin a — Suf^x cos 

011 = ^ a — 5 Tify X cos (T^ 
b 

/3i2 = — ^ ( —5 X cos cr + 6 rtfy x sin a J 

/3i3 — ^ + (5 X sin cr^ 



/3j^4 = ■g-^ (""^ ^ o" — <5 nj, X cos cr^ 



Vi(t) = 1 

l/'2(t) = " 
V3(t) = «^ 

V'4(t) = 

*5(t) = 

^6(t) = 

ljjj(t) — COS u 

ipS (i) — sin u 
ijjg (t) = u sin u 

■010 (*) — WCOS'U, 

V'li (*) = cos w 

''/'12 (*) = sin u 

V'lsC*) = it^sinu 

V'14(*} = cos W 



are significant in fitting the rotational and orbital phases of the pulsar over the available time span 

of observations. 

Let us introduce auxilary functions defined according to the formula 



14 



6=1 

where the matrix of information 



(19) 



niN 



(20) 



i=l 



the matrix L^^ is its inverse, and AT = NPi, is a total span of observational time. Functions E{t) 
are dual (Deeter 1984) to ■0a(^) because of the cross-orthonormahty condition 



mN 



J2'^a{ti)iJb{ti) = Sab- 



(21) 



i=l 
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Now suppose that we measure m equally spaced and comparably accurate arrival times for each 
orbit for a total of N orbital revolutions, so we have mN residuals = r(tj), i = l,...,mN. 
Standard least squares procedure (Bard 1974) gives the best fitting solution for the estimates of 
the parameters Pa 



mN 



Pa{AT) = Y,Ea{ti)e{ti), a = l,...,14. (22) 

1=1 

Let the angular brackets denote an ensemble average over many stochastic realizations of the 
observational procedure. Hereafter, we assume that the ensemble average of the noise e{t) is equal 
to zero. Hence, the mean value of any parameter (3a is equal to zero as well, i. e. 

< e(t) >= 0, and < (3a >= 0. (23) 
The covariance matrix Mat = < (3a(3b > of the parameter estimates is given by the expression 

mNmN 

Mab{AT) = J2J2^a{U)E,{tj)R{t^,t,), (24) 
i=lj=l 

where R{ti,tj) = < e(ti)e(tj) > is the autocovariance function of the stochastic process e(t). The 
covariance matrix is symmetric (that is, Mab = Mba), elements of its main diagonal give variations 
(dispersions) of measured parameters cr^^ = Maa = < Pi >, and the off-diagonal terms represent 
the degree of statistic covariance (correlation) between them. The covariance matrix consists of two 
additive components and describing contributions from the non-stationary and stationary 
parts of autocovariance function R{ti,tj) respectively. Explicit expressions for the matrix Mab can 
be found in the paper (Kopeikin 1999b) wherein we have done all our calculations in the time 
domain. Only admits the transformation to the frequency domain and this Fourier transform 
will be discussed in subsequent sections along with particular spectral properties of M^. 

Subtraction of the adopted model from the observational data leads to the residuals which 
are dominated by the random fluctuations only. An expression for the mean-square residuals after 
subtracting the best-fitting solution for the estimates ()22|) is given by the formula 

1 mNmN 

< r\AT) >= —^T.T.K{U,t,)R{t,,t,), (25) 

"'^^ i=ij=i 

where the function 

14 

K{t,, tj) = 5i, - Y.'Ea{ti)iJa{tj), (26) 
a=l 

is called the filter function (Blandford et al. 1984). We have proved (Kopeikin 1999b) that the 
post-fit residuals depend only on the stationary part of the noise (even in the case when the non- 
stationary part of the autocovariance function is present) which means that Eq. (|25|) can be written 
down as 



a=lfe=l 



14 14 

mN : 



mNmN 

J2Y.Mu)Mti)R'iu,t, 
i=ij=i 



(27) 



For this reason, methods of spectral analysis in frequency domain can be applied for analyzing 
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residuals without any restrictions. Let us note that the exphcit dependence of TOA residuals on 
the total span of observations is contained in (Kopeikin 1999b). 

4 FOURIER TRANSFORM OF FITTING FUNCTIONS 

We define the Fourier transform of the fitting functions ipait) as 

mN 

^aif, m, N) = expi-2mft,), (28) 

i=i 

where / is the Fourier frequency measured in units inversly proportional to the units of measurement 
of time t. We measure time in units of the orbital phase u = ribt, that is in radians. Then the 
frequency u = 2-7?/ is dimensionless and is measured in units of the orbital frequency n?,. One notes 
the Fourier transform of the fitting functions depends on three arguments: the Fourier frequency 
/, the total number of orbital revolutions A^, and the rate of observations m. 

When the total number of observational points, mN, is large enough we can approximate the 
sum in Eq. (|28j) by the integral (Kopeikin 1999b) 

~ m ~ 

^a(c^,m,Ar) = — V'a(cu,Ar), (29) 
'ijja{uj,N)= / ijjaiu) exp{—iu!u)du. (30) 

J-ttN 

We note that ipa{—^^) = i^li^^), where the asterisk denotes a complex conjugation. Replacing the 
sum over observational points by the integral with respect to time (the orbital phase) makes our 
calculations equivalent to the case of continuous observations with uniform distribution over time. 
The following formulae are also of use in practical computations 



rnN 

2 ^paiu) cos{uju)du, if index a = 1, 3, 5, ... 
Jo 

(31) 

pttN 

—2i ipaiu) sm(uju)du, if index a = 2, 4, 6, .... 
Jo 



These expressions shows that the fitting functions with odd indices are real and those with even 
ones are complex. 

Let us denote T = vrA^ and z = uT. One notices that T relates to the total span of observation 
using equation T = nhAT/2. The Fourier transform of fitting functions takes the form 

^i(^) = 2T0i(^), (32) 

i>2ito) = 2tT^2iz), (33) 

M^) = 2T303(^), (34) 

ip^itu) = 2iT^04(2), (35) 

^,{uj) = 2T'Mz), (36) 

i^eiuj) = 2«T606(^), (37) 
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i/;q{uj) 
-011(0; 

-012 (CU 

-014(0; 



^i(^ + T)+0i(^-T) 



iT^ 



iT^ 



03(^ + T)-03(^-T) 

04(^-T)+04(-2 + T) 



(38) 
(39) 
(40) 
(41) 
(42) 
(43) 
(44) 
(45) 



where functions (pai^) (a = 1, 2, 6) have the following form 



(46) 



cos z sm z 



(47) 



sin z 2 cos z 2 sin z 



(48) 



7 cosz Ssinz 6cosz 6 sin 2; 

M^) = — -o 1^ + ^^. 



(49) 



05 (^) 



06 (^) 



sin z 4 cos z 
+ ^ 



12 sin^ 



Z'' 



24 cos z 24 sin z 
^ + 



cos 2; 5 sin 2; 20 cos 2; 60 sin 2; 120 cos 2; 120 sin 2; 
h , \ z . 



(50) 



(51) 



z-" z-^ z^ 

Simple inspection reveals that these functions are linear combinations of the spherical Bessel 
functions jjyz) defined as (Korn & Korn 1968, section 21.8-8) 

Id'" 



sm z 



z dz 



(a = 0,1, 2,...). 



(52) 



The Bessel functions have a polynomial behaviour ^ near the point 2; = 0, and then oscillate with 
monotonically decreasing amplitude. Asymptotic expansion of the sperical Bessel functions for large 
values of the variable z (that is directly proportional to the Fourier frequency) are given by the 
formula 



4 The function ja{z) = + 0(2;"+2) for z < 1 (a = 0, 1, 2, ...). 
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sm I z 

Ja{z) ^ \ ' ■ (53) 

It is worth emphasizing that the maximum value for any of these functions can not be larger than 
1. 

The fitting functions (f)a{z), expressed in terms of the spherical Bessel functions, assume the 
form 



(pi{z 



Uz), (54) 

-Ji{z), (55) 
1 2 

^joiz) - ^j2{z), (56) 

-Ijiiz) + lUz), (57) 

lUz) - ^Mz) + ^J4{z), (58) 

-^Jiiz) + ^Mz) - ^^3,{z). (59) 



5 FOURIER TRANSFORM OF THE COVARIANCE MATRIX 

In order to calculate the covariance matrix we need to know the Fourier transform of the dual 
functions Sa(t). The transform is defined in accordance with definition (fT^ of the dual functions 
and takes the form 

S„(/, m, AT) = Y.Kl^c{f. m, AT), (60) 

c=l 

and the cross-ort honor mal condition in the frequency domain is given by the integral 

H,(/, m, iV)§,(/, m, N)df = ^5,6. (61) 
2 

In the limit of continuous observations it is convenient to introduce the matrix Cab = ^^ab 
instead of the information matrix Lab- The explicit expression for the matrix Cab is given by the 
integral: 

/nN 
ilJa{u)ipb{u)du, (62) 
-nN 

and the result of evaluation of this integral is given in the paper by Kopeikin (1999b) in Tables 5 
and 6. Then we have L~^^ = (27r/m)C^^, and the Fourier transform of the dual function S(/) can 
be recast as 

^a{f,N) = j:C-,'Mf,N). (63) 

b=l 
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Hence, comparing the eq. ()63|) with ()60|) one concludes that in the hmit of continuous observations 
the Fourier transform of the dual functions depends only on the Fourier frequency and the total 
number of orbital revolutions as was expected. It is more insightful to express the dual functions 
(jF)!^ in terms of the spherical Bessel functions (jK^ . The expressions obtained are rather unwieldy, 
and, for this reason they are given in Appendix A. 

Making use of the definition of the Fourier transforms of the stationary part of the autocovari- 
ance function (jSJ and the dual functions (jUnj) we obtain the stationary part of the covariance matrix 
M~i^{m,N) (see Eq. (j2D) for its definition) expressed as follows 



M^, = J^ S{f)HM,m,N)df, (64) 
where Hab{f,iTL,N) is the transfer function given by the expression 

Habif, m, N) = Eaif, m, N)El{f, m, N) + E^f, m, N)E,{f, m, N), (65) 

and an asterisk denotes complex conjugation. For numerical computations of the following 
formula can be used in practical computations 

M-,{m,N) = ^ f\^,{f,m,N)f--df+ (66) 



(2vr) 

Ik 



e 



where derivatives of Hab are taken with respect to the Fourier frequency, ellipses denote terms of 
negligible influence on the result of the computation, and A is the upper cutoff frequency arising 
from the sampling theorem and is inversly proportional to the minimal time between subsequent 
observational sessions. We emphasize that the matrix M~f^{m,N) does not depend on the lower 
cutoff frequency because all contributions from the lower limit of the integral in Eq. ()66|) are 
canceled out by corresponding terms from the series. 



6 FOURIER TRANSFORM OF THE FILTER FUNCTION 



Timing residuals are expressed in terms of the Fourier transform of the filter function. Employing 
eqs. (0), (j2ZI), and (I^Uj) yields for the residuals 



oo 



<r'>=2 S{f)K{f,m,N)df, (67) 



where is the Fourier transform of the filter function (|26|) 

Kif, m,N) = l- [h.(/, m, N)^:if, m, N) + E^f, m, N)^M, ^, N)] , (68) 

is the Fourier transform of the filter function defined in Eq. (jSEI)- In the limit of continuous obser- 
vations there is no dependence on the frequency of observations, m, so that one obtains 

K{f, AT) = 1 - i-^ [em, N)ra{f, N) + Elif, N)Mf, N)] . (69) 

Plots of the Fourier transform of the filter function K{f) are shown in Appendix D for 
different numbers of orbital revolutions A^. In any case the filter function is approximately equal to 
1 until the frequency is higher than 1/(AT), and rapidly decreases in amplitude as the frequency 
approaches zero. We also note that the plot of the Fourier transform of the filter function clearly 
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shows the additional dip near the orbital frequency. The dips near zero and orbital frequencies get 
narrower as the number of observational points increases. 



7 SPECTRAL SENSITIVITY OF TIMING OBSERVATIONS OF MILLISECOND AND BINARY 
PULSARS 

Analytical expressions and graphical representations of Fourier transforms of dual functions and 
timing residuals help us to understand in more detail the spectral sensitivity of single and binary 
pulsars to different frequency bands in the spectral decomposition of the noise. First, let us consider 
behaviour of the Fourier transform of fitting functions near zero and the orbital frequencies. 

It is easy to confirm after making use of Taylor expansion of exponential function in near 
~ that 



if a=l,3,5,... 



Cal — 2^ Ca3 + + ^ Pa, 

i (^-ujCa2 + luj^Cai - i^^^C^e + uj'^Pa) , ^ a=2,4,6, 



(70) 



where pa is a residual term depending only on the total number of orbital revolutions, N. A Taylor 
expansion of the Fourier transform of the fitting functions near the orbital frequency yields 



if a=l,3, 



-Ca8 +{UJ- l)a.lO + 



Ca.l2 — 



Ca.ii + (t^ - lYqa , if a=2,4,..., 



(71) 



where qa is a residual term depending only on the total number of orbital revolutions, (recall 
that the frequency is measured in units of the orbital frequency rih). Applying to Eqs. (fTU jl -ffTT I) 
the definition of the dual functions ()63|) in the limit of continuos observations yields the asymptotic 
behaviour of the dual functions 



if a=l,3,5,... 



i (-u;5a2 + |tu^5a4 " j^^^^ad + tu'^^a) , if a=2,4,6,.... 



(72) 



near zero frequency, and 



5a.l3 + iuj - l^Qa, 



if a=l,3, 



-6aS + (^ - l)Sa.lO + ^^5a.l2 " ^^5a.l4 + (c^ " l^Qa] , if a=2,4,.. 



(73) 



near the orbital frequency, where Pa and Qa are residual terms. Table 2 shows the asymptotic 
behaviour of the residual terms of the dual functions. 

Now we can study the asymptotic behaviour of the filter function K{f) as defined by eq. 
Taking into account the fact that Y^lti CatCbc = Sac where 6ac is the unit matrix, we obtain 



K{f) 



bo ■ uj^ 



when io 



bi ■ (cu — 1)^, when u ^ 1 



(74) 



where bo and bi are numerical constants which can be calculated precisely but are not important 
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for the discussion present in this section. Such a dependence of the filter function K{f) on the 
frequency / significantly reduces the amount of the detected red-noise power below the cutoff 
frequency /cutoff — ctcT^^ and in the frequency band 1 — a^T^^ ^ / ^ 1 + a^T^^ lying near 
the orbital frequency. Here constant coefficents ac and ab can be determined by comparing the 
calculations of the mean value of the timing residuals in time and frequency domains (Kopeikin 
1997a). The low Fourier frequencies in the noise power spectrum are fitted away by the polynomial 
fit for the spin-down parameters of the observed pulsar. The Fourier frequencies close to the orbital 
frequency are fitted away by the fit for the orbital parameters of the pulsar. The amount of noise 
power remaining in timing residuals after completion of fitting procedure ^ is estimated by the 
expressions 

<r^>=2/^ S{f)df + 2 S{f)df 
and 

<r' >=2r^S{f)df + 2r S{f)df 

where for the calculation of the first and the second integrals we have taken the very first terms 
of the spectra ^ given in Eqs. (|T0|l - (fT5|l and assumed that a^/T ^ 1 and a^/T ^ 1 so that these 
quantities are used as small parameters of the Taylor expansion of the integrals under calculation. 
We draw attention of the reader that in real signal processing the infinite-frequency limit in Eqs. (j75|l . 
()7(ij) is, in fact, inversely proportional to the sampling time of observations. 

The second term in the right hand side of Eqs. (f73j) . ()76|) shows amount of noise absorbed by 
fitting orbital parameters. It is negligibly small compared to the first term in the right hand side 
and can be neglected in practice. Hence, we conclude that the post-fit timing residuals can be used 
for the estimation of the amount of red noise and its spectrum in frequency band just from ttcT"^ 
up to infinity, irrespectively of whether the pulsar is binary or not. This reasoning puts on firm 
ground the estimates of spectral sensitivity of timing observations and the cosmological parameter 
fig, characterizing energy density of stochastic gravitational waves in early universe, made by Kaspi 
et al. (1994) and Camilo et al. (1994) using observations of binary pulsars PSR B1855+09 and PSR 
J1713+0747 respectively. 

Analysis of the spectral sensitivity of the estimates of variances of spin-down and orbital pa- 
rameters is more cumbersome. We are interested in which frequencies give the biggest contribution 
to the variances. This is important to know, for example, if we want to use variances of certain 
orbital parameters for setting the fundamental upper limit on Qg in the ultra-low frequency domain 
(Kopeikin 1997a). Analitical calculations reveal the leading terms in the asymptotic expansions of 
the dual functions near zero frequencies which are present in Table 2. The squares of the dual func- 
tions appearing in eqs. ()B1|) - ()B14|) and eqs. ()C1|) - ()C14|) are rather complicated. Their behaviour is 
periodic with bumps both near zero and the orbital frequencies and with occilating behaviour which 
amplitude is rapidly decaying far outside of these frequencies. We determined that the amount of 
noise absorbed by fitting of the spin-down parameters near zero frequency is essentially bigger than 
that near the orbital one. On the other hand, the amount of noise absorbed by fitting of the orbital 
parameters is substantial for the Fourier freqiencies lying near the orbital frequency. Thus, there 
are two spectral windows in which parameter' variances of the timing model are most sensitive to 
the stochastic red noise and they are located near zero and the orbital frequency. These windows 




^ It is useful to compare calculations of the residuals given in the present paper with those given in (Kopeikin 1999b. 
® Delta-functions and their derivatives do not contribute to the integrals in the approximation under consideration. 
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Table 2. Asymptotic behaviour of residual terms of the dual functions near zero and orbital 
frequencies. Constant h = cosT = (—1)^. 



Dual function Residual term Pa Residual term Qa 



rp6 l^T^h 

r ^ 56 



-iT^h 



1584 4 



6864 4 

l_rp2 45 h 

528 8 

1 rp2 33 h 

3120 40 

1 q^4u 1 a-'4 

165 280 



_rn6u ^X'' 

5 280 



693 28 



273 28 



i ^ " 28 



^—T^h 5-T2 

9009 28 

^T^h i 

297 12 



3861 12 



are restricted by two frequency intervals, (0,^), and, (1 — 1 + ^), respectively, where constant 
coefficients a, a_, a+ can be calculated by comparing calculations of variances of the fitting pa- 
rameters in time and frequency domains. The variances of the fitting parameters depend only on 
the total span, AT, of observations (Kopeikin 1997b). Comparing the dependence of the variances 
on AT calculated in time domain with that calculated in the frequency domain with the integrals 
truncated by the two frequency windows, reveals the relative importance in contribution of different 
frequencies to the integrated values of the variances of the parameters given by two integrals 



Ii ~ hr 



and 



2k 



df, 



(77) 



I9 



(2 



1+- 



T 1^ 



(7J 



The ffist integral describes the contribution of low frequencies to the parameter's variances. 
The second integral gives the contribution to the parameter's variances from the frequencies lying 
near the orbital frequency. It is worth emphasizing that the terms with delta functions and their 
derivatives in Eq. ()77|) cancel out all terms depending on the cutoff frequency e and diverging as 
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Table 3. Comparative contribution of different frequency bands to variances of spin-down (a = 
1,2, ...,6) and orbital (a = 7,8, ...,14) parameters I3a- Number n = 1,2, ...,6 denotes the spectral 
index of corresponding red noise. Time dependence of all variances completely coincides with that 
which was obtained by calculations in time domain as given in (Kopeikin 1999b). 



Variance of parameter Contribution of integral Ii Contribution of integral I2 



„2 


^ — 1 


^ T"-3 
rsj 1 


J2 


~ 1 


, , rp — 5 




^ , T-'n — 5 




„2 


rr<n — 7 
~ 1 


rp-9 




^ rpn-9 


^ T-11 


'^k 




~ T-13 


'^k 




^ T-i 


'^k 




^ T-i 


'^k 




^ T~3 






^T-3 


pii 




^ T~5 






^ T~5 




^ '-£71 — 9 


^T-7 


"k. 


^ rj-in — 11 


~T-7 



£ goes to zero. This cancellation has been expected since we modified the spectrum of red noise 
to avoid the appearance of all divergent terms which have no physical meaning. We don't give 
here the results of the calculation of numerical values of constants a, a_, a+ because they are not 
so important for general conclusions. The asymptotic behaviour of the dual functions near zero 
and the orbital frequency is enough informative to see which frequncy band is the most important 
for giving contributions to the corresponding integrals and the parameter's variances. The time 
dependence of two integrals is shown in Table El up to a not so important constant. The behavior 
of the variances of the first three spin-down parameters is not physically interesting because they 
are contaminated and strongly biased by the presence of the non-stationary part of the red noise 
(Kopeikin 1997b). For the rest of the spin-down parameters one observes that the contribution of 
the noise energy from low frequencies to the variances of the parameters is dominating. However, 
the situation is not so simple in the case of the orbital parameters. One can see that in the case 
of red noise having spectral index n ^ 4, the contribution of the noise energy from the orbital 
frequency interval (1 — ^,1-1- ^) can be equal to or even bigger than that from the low frequency 
band. Only when the spectral index of noise is n ^ 5 does the contribution of the noise energy of 
low frequencies to the variances of orbital parameters begins to dominate. 

It is worth noting that the timing noise with the spectral index n = 5 is produced by the 
cosmological gravitational wave background. The fact that for this noise low-frequencies give the 
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main contribution to the variances of orbital parameters confirms our early statement (Kopeikin 
1997a) that the measurement of the variances of orbital parameters showing secular evolution 
probes the ultra-low frequency band of the cosmological gravitational wave background. Hence, 
these variances can be used for setting an upper limit on the cosmological parameter Qg in this 
frequency range, in contrast to timing residuals which test only the low-frequency band of the 
background noise as explained in (Kopeikin 1997a). 
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APPENDIX A: EXPLICIT EXPRESSIONS FOR THE DUAL FUNCTIONS 

In this appendix we give explicit expressions for the dual functions. Using definition ()63|) and 
elements of inverse matrix from the paper (Kopeikin 1999b, Tables 5 and 6) one obtains 



5 27 

Si (2) = Mz) + -j2{z) + yj4(2) + 



-^{3 [Ji{z + T) - j,{z - T)] - 7 [Uz + T) - Uz - T)]} 
-^{3 lUz + T) +Uz - T)] - 20 [Uz + T) +Uz - T)]} 



(Al) 



-22W 



3 
T 



ji(^) + ^Mz) + ^j^iz) 



+ 



^|jo(^ + T) - jo(^ - T) - 5 [j,{z + T)-Uz-T)]j + 
^jsi [Mz + T) + j,{z - T)] - 154 [Mz + T) + Mz - T)]| 



(A2) 



23(^) 



15 



2T2 
105h r 



9 

Mz) + -jiiz) 



1 ()5h r ^ 

^{3 [Jliz + T) - n{z - T)] - 7 [Uz + T) - Uz - T)]| 

in5h r ^ 

[jo{z + T) + Uz - T)] - 50 [Uz + T) + U^ - T)]| , 



(A3) 



35 
2T3 



Jsiz) + —Uz) 



(A4) 
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-^{M^ + T) - Mz - T) - 5 \j,{z + T) - Mz - T)] } 

' |9 \j^{z + T) + Mz - T)] - 28 \j,{z + T) + h{z - T)] | , 



4X5 



+ ^{3 [Ji{z + T) - - T)] - 7 [Mz + T) - j3(^ - T)]| 
— ^{bo(;^ + T) + Mz - T)] - 8 b2(;. + T) + M^ " T)]| , 



1^ / N 693 ., , 



693h f 
"8T6 

2079h r 



|jo(^ + T) - jo(^ - T) - 5 \j2{z + T) - - T)]| + 
^{13 Mz + T) + - T)] - 42 b3(^ + T) + M^ - T)]} , 
S7(^) = io(-2 + T)+io(-2-T) + ^[72(^ + T)+j2(-2-T)] 

-^{^ b-i(^ + T) - - T)] - 7 [Mz + T) - j3(^ - T)]} 

b2(^) - 6j4(^)] , 



^2 



-Es{z) = M^ + T) - M^ - T) + ^ b2(^ + T) - j2(^ - T)] + 

+ T) + Mz - T)] - 7 [j^iz + T) + j3(^ - T)]} + 
Y[3JiW-7J3W + 11J5W] , 
S9(^) - |{ji(^ + T)-ii(z-T) + ^L73(^ + T)-i3(^-T)]} 

-{ [jo(^ + T) + io(^ - T)] - 5 L72(^ + T) + Mz - T)] } 
[Mz) - 5Mz) + 9Mz)] , 



15 
1X2 

15h 



^2 



-fioiz) = -|{bi(^ + T)+ji(;.-T)] + ^b3(;^ + T)+j3(^-T)]} 
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15 



4T2 
15h 



|io(^ + T) - jo{z - T) - 5 \j2{z + T) - j2{z - T)] + 



^3 



[-18ji{z) + 77n{z)- 220j,{z)] , 



Sii(^) = -^,Mz + T)+j2{z-n (All) 



15 

+ 



4T3 
15h 



|3 \j,{z + T) - j,{z - T)] - 7 L73(^ + T) - Mz - T)]} + 



'J'4 



[2io(-2)+5i2(^)-72j4(^)] , 



7Si2(^) = -^b2(;. + T)-j2(;.-T)] (A12) 



15 



4T3 
15h 



{3 + T) + j,{z - T)] - 7 b3(^ + T) + Mz - T)] I + 
[3ji(^)-7j3(^) + llj5(^)] , 



35 

Si3(^) = -^b3(^ + T)-i3(^-T)]+ (A13) 



35 
4T4 

35h 
^4 



{jo(^ + T) + joiz - T) - 5 b2(^ + T) + j2(^ - T)] I + 

IMz) - 5j2iz) + 9u{z)] , 



1 35 

t5i4(^) = ^b3(^ + T)+j3(^-T)]+ (A14) 



35 
4T4 

105h 



{jo(^ + T) - jo(^ - T) - 5 \j2{z + T) - j2(^ - T)] I + 
[4ji(^)-21j3(^)+66j5(^)] . 



APPENDIX B: ASYMPTOTIC BEHAVIOUR OF THE DUAL FUNCTIONS NEAR ZERO 
FREQUENCY 

In this appendix we give the asymptotic behaviour of the dual functions near zero frequency. They 
are as follows: 

~ ~ ~ 5 27 

= Uz), ^,(z) ^ jo{z) + -Uz) + -u{z) , (Bl) 
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1 ~ ~ 7 'I'l 

^S,(/) = --Uz), U^) ^ ji{z) + ^js{z) + -Mz) , (B2) 

S3(/) = Uz)^j2{z) + lMz), (B3) 

1h4(/) = ^e4(^), Uz)-js{z)+^-^Mz), (B4) 

55(/) = e5(^)=j4(^), (B5) 

= -^M')^ Uz)-h{z), (B6) 

= -^Uz), ir{z)^j2{z)-6j,iz)-^zsmz, (B7) 

1~ 9h~ ~ 7 11 1 

tH8(/) = Y^'^^^' ^8(^)=ii(^)-3i3(^) + yi5(^)-3smz, (B8) 

15h ~ 

E,{f) = -^Uz), Uz)^jo{z)-5j2iz) + 9u{z)-cosz, (B9) 

1~ 270hr , , ~ ^ , , 77 , , 110 , , 5 1 , 

-^io(/) = — ^6oU) , Cio(^) = Ji(^) - Y^Jsl^) + -^J5(^) - Y^sm^- — ^cos^ , (BIO) 

30h ~ ~ 5 1 

Sii(/) = ^^11(2;) , ^11(2:) = Mz) + ^j2{z) - 36ji{z) - -zsinz- cos 2; , (Bll) 

]^i2{f) = -^Uz), Uz)=Mz)-lMz) + ^Mz)-^smz, (B12) 

Si3(/) = i^63(;^), 63(^) = 6(;^), (B13) 

1~ 420h~ , , ~ ^ , , 63 , , 33 , , 1 1 

T^i4(/) = ^^Ci4(^) , Ci4(^) = Ji{z) - Y^j3{z) + yjsl^) - -sm^ - —zcosz , (B14) 

where h — (—1)-^. 

APPENDIX C: ASYMPTOTIC BEHAVIOUR OF THE DUAL FUNCTIONS NEAR ORBITAL 
FREQUENCY 



The leading terms in the asymptotic expansions of the dual functions '^a{z) near the orbital fre- 
quency are as follows: 

45h~ ~ 7 1 

= --^6(2/), Ci{y) = Ji{y) - ^jsiy) - ^ sin y , (CI) 

-^2(/) = -^Uv), 6(y)=Jo(2/)-5j2(2/)-cos|/, (C2) 
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53(/) = -^e3(y), i3iy)^Uy) (C3) 

]Uf) = -yMv)^ Uy) = Uy), (C4) 

55(/) = §e5(z/), e5(?/)=6(?/), (C5) 

i^eU) = ^e6(z/), Uy)-Uy), (C6) 

57(/) = My), i7{y)^jo{y) + lj2{y), (C7) 

^S8(/) = e8(z/), Uy)^-i7{y), (C8) 

59(/) = -|e9(y), iM^JM + ljsiy), (C9) 

-5io(/) = ;^eio(y), 6o(y)=6(2/), (CIO) 

5ii(/) = -^eii(y), 6i(y)=j2(y), (cii) 

^5i2(/) = ;^ei2(y), Ii2(y) = -lii(y) , (Cl2) 

35 

Si3(/) = ^ei3(y), 63(y)=j3(y), (Ci3) 

tSi4(/) = ;^li4(z/), Ii4(y) = Ii3(y) . (C14) 
where y — z — T. 

APPENDIX D: PLOTS OF THE FOURIER TRANSFORM OF AUTOCOVARIANCE FUNCTION 
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Figure Dl. Plot of the Fourier transform of tlie filter function of timing residuals for the amount of orbital revolutions N = A. Frequency 
is measured in units of orbital frequency nj . Amplitude of the transform has been normalized to unity. 




Figure D2. Plot of the Fourier transform of the filter function of timing residuals for the amount of orbital revolutions N = 10. 
Frequency is measured in units of orbital frequency n;,. Amplitude of the transform has been normalized to unity. 
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Figure D3. Plot of the Fourier transform of tire filter function of timing residuals for the amount of orbital revolutions N = 30. 
Frequency is measured in units of orbital frequency Uf,. Amplitude of the transform has been normalized to unity. 



